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ABSTRACT 


Kozai-Lidov (KL) oscillations in hierarchical triple systems have found application 
to many astrophysical contexts, including planet formation, type la supernovae, and 
supermassive black hole dynamics. The period of these oscillations is known at the 
order-of-magnitude level, but dependencies on the initial mutual inclination or inner 
eccentricity are not typically included. In this work I calculate the period of KL oscil¬ 
lations (£rl) exactly in the test particle limit at quadrupole order (TPQ). I explore 
the parameter space of all hierarchical triples at TPQ and show that except for triples 
on the boundary between libration and rotation, the period of KL oscillations does 
not vary by more than a factor of a few. The exact period may be approximated to 
better than 2 per cent for triples with mutual inclinations between 60° and 120° and 
initial eccentricities less than ~0.3. In addition, I derive an analytic expression for 
the period of octupole-order oscillations due to the ‘eccentric KL mechanism’ (EKM). 
I show that the timescale for EKM oscillations is proportional to e oct , where e oc t 
measures the strength of octupole perturbations relative to quadrupole perturbations. 


1 INTRODUCTION 

Triple systems are common in the Galaxy, compris¬ 
ing ~10% of systems in whic h the primary i s ~1 
Mm dDuauennov fe Mavoi _ 1991; iRaghavan et al.1 l20ld : 


iTokovinin 20141 : Riddle et al.l 201f ). All observed triples are 
‘hierarchical,’ in that the relative distance between two com¬ 
ponents of the triple is much smaller than the relative dis¬ 
tance between them and the third 1 S uch systems are sta ble if 
they are sufficiently hierarchical dMardling fe Aarsethll 19991 . 

booifi FI 

In general, if the tertiary is highly inclined relative to 
the inner binary, the eccentricity of the inner binary will un¬ 
dergo oscillations, known as Kozai-Lidov (KL) oscillations 
jLidovlll962l :lK oza KL oscillations have been invoked 

in many contexts to explain a wide v ariety of phenomena 
such as the formation of hot Jupiters dWu fe Murravll 2003: 


20071: Fabrvck v fe Tremainel 20071 : Naoz et al.1 


Wu et al.l _ ... _ _ _ __ _ _ 

201ll . 12012: Petrovichl 20151). the formation of blue strag¬ 


glers dPerets fe Fabrvckvll2009l: |Naoz fe Fabrvck vll2014l). the 


merg er of WD-WD binaries dThompsonl 201ll : iKatz et al.l 
120111) . the merger of supermassive and intermediate- 
mass black holes dMiller fe Hamiltonl l2002bl : iBlaes et al.l 


1 While sta ble non- h ierarchical trip le systems are possible (e.g ., 

IChenciner fe Montgomervll2000l ; rSuvakov fc D mitraginovidfipiJ) . 

they may require fine tuning to form and have never been ob¬ 
served in nature. 


120021 : I Wenl [20031 ). the distribution of dark matter around 
supermassive black hole binaries dNaoz fe Silkl 120141 ) . 
and as a source of unique gravitational wav e sig¬ 
nals dM iller fe Hamiltonl 2002a: Goul dl 20 111: Setd 20131: 


Antonini et al.l 120141 : lAntognini et al.l 120141 : iBode fc Wega 


20141) . 


KL oscillations are a secular phenomenon, occurring on 
timescales much longer than the orbital periods. It is there¬ 
fore possible to average the motions of the individual stars 
over their orbits and study only the secular changes to the 
orbital elements. If there is a large mass ratio in the inner 
binary then on even longer timescales the strength of the 
KL os cillations (i . e., th e maximum eccent r icity reached) will 
vary dFord_et_jilJ 2000l : IKatz et al.l 1201 ll : iLithwick fe Nao3 
l201ll : lNaoz et af1l2013al ). These variations have been termed 
the ‘eccentric KL mechanism’ (EKM), and in some cases can 
cause the inner binary to pass through an inclination of 90° 
with respect to the outer binary in a ‘flip’ from prograde 
to retrograde or vice versa. During a flip the eccentricity of 
the inner binary can be driven to extremely large values be¬ 
cause the strength of KL oscillations is very sensitive to the 
mutual inclination, with arbitrarily strong oscillations oc¬ 
curring as the inclination approaches 90° exactly in the test 
particle limit. Although EKM oscillations do not occur when 
the two stars of the inner binary are of equal mass, mass 
loss from one of the stars in t he course of st ell ar evolution 
can induce EKM oscilla tions dShappee fe Thompsonl [20131 : 
iMichaelv fe Perefisl 1 201 4l) . EKM oscillations and flips have 
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2 Antognini 


generally been studied in the context of hierarchical triples, 
but flips occur over a wider range of p arameter space in 
both 2+2 quadruples (IPeicha et al.ll2013l f and 3+1 quadru¬ 
ples llHamers et al.ll2015l b 

Because the extreme eccentricity oscillations that occur 
during a flip can affect the evolution of the objects in the 
inner binary, the timescale for EKM oscillations is another 
important quantity. Yet no derivation of the timescale for 
EKM oscillations has appeared in the literature, although 
several studies have as serted that Iekm ~ fKL/eoct is a 
plausible timescale fe.g-. lKatz et al.ll201ll ; lNaoz et al.ll2013bl : 
Li ct al .1 120150 . where e oc t measures the strength of the oc- 
tupole order term relative to the quadrupole order term of 
the Hamiltonian (see equation 1501 for a definition). I show 
that Iekm ~ Ikl / y/Coct- 

In Section [2] I present the basic parameters and equa¬ 
tions that govern a hierarchical three-body system. In Sec¬ 
tion [3] I then derive the period of KL oscillations. In Sec¬ 
tion [4j I explore how the period varies over the parameter 
space and in Section [5] I provide an approximation to the 
exact period. In Section [U] I treat the period of EKM os¬ 
cillations and derive the corrected timescale. I conclude in 
Section [7] 

To perform the calculations in this paper I wrote the 
kozai Python module. This module can evolve hierarchical 
triple systems in the secular approximation up to hexade- 
capole order using either the Delaunay orbital elements or 
the eccentricity and angular momentum vectors. I have re¬ 
leased this code under the MIT license and it is available at 
https://github.com/joe-antognini/kozai. 


2 BASIC EQUATIONS 
2.1 Notation 

Throughout this paper orbital properties referring to the 
inner and outer binary are labelled with a subscript 1 and 
2, respectively. The masses of the two components of the 
inner binary are mi and m 2 , and the mass of the tertiary is 
m 3 . 

We will often refer to the orbital parameters using De¬ 
launay’s elements: the mean anomalies, l x \ the arguments 
of periapsis, g x , and the longitudes of ascending nodes, h x , 
where x — 1 or 2 and refers to the inner or outer binary, 
respectively. Their conjugate momenta are 

Li = mim2 G(mi + m 2 )a!, (1) 

mi + m 2 


L 2 


m3 (mi + m 2 ) 
mi + m 2 + m3 


V G(mi + m 2 + m 3 )a 2 , 


( 2 ) 


G x = L X \J 1 — e 2 , (3) 

and 

H X — G X cos i x * (4) 

Delaunay’s elements form a set of canonical variables. Note 
that G 1 and G 2 are the angular momenta of the inner and 
outer binaries, respectively. We furthermore define the re¬ 
duced angular momentum 

jl = l-el. (5) 


The angular momentum of an orbit may thus be written 
G x — L x j x . 


2.2 The Hamiltonian 


If a three-body system is sufficiently hierarchical, its Hamil¬ 
tonian may be considered to be that of two isolated binaries 
(the inner binary, consisting of the two closest bodies, and 
the outer binary, consisting of the distant body plus the 
inner binary taken as a point mass) plus a perturbative in¬ 
teraction term: 


_ Gmim 2 G(mi + m 2 )m 3 

iL — ^ 1 1 /Tpert- 


2oi 


2a 2 


( 6 ) 


This interaction term captures the change in the orbital mo¬ 
tion of each binary in the tidal field of the other. Because we 
are assuming that the triple is hierarchical, the semi-major 
axis ratio, a — ai/a 2 , is a small parameter that we can use 
to expand the perturbati ve component of the Hamiltonian 
in a multipole expansion (IHarringtonll 19681 1. 


T^pert 


„ 00 




Pj{ cos4>), 


(7) 


where Pj is the j th Legendre polynomial, r x is the distance 
between the two components of the binary, $ is the angle 
between r 2 and n and Mj is a mass parameter defined by 


+-1 


K A 

Aij = 17111712 ^ 3 - 


- (~m 2 ) 


3 — 1 


(mi + m 2 y 


( 8 ) 


If we are only interested in changes to the orbital el¬ 
ements that occur on timescales much longer than the or¬ 
bital periods (so-called ‘secular’ changes), we must average 
the Hamiltonian over both mean anomalies. To do this while 
maintaining the canonical structure of the Hamiltonian re¬ 
quires a technique known as von Zeipel averaging. The gen¬ 
eral case for three massive bodies is quite complicated even 
at quadrupole order as one m ust b e careful to include the 
longitudes of ascending nodes llNaoz et al.ll2013ah . However, 
the Hamiltonian simplifies considerably if one component of 
the inner binary is taken to be a test particle as this allows 
one to fix the longitudes of ascending nodes and eliminate 
them from the Hamiltonian. The resulting double-averaged 
Hamiltonian at quadrupole order in the test particle limit is 


Hq = C 2 [(2 + 3e?) (l — 3 cos 2 i) — 15e? (l — cos 2 i) cos 2gi] , 

(9) 

where G 2 is a constant parameterizing the strength of the 
quadrupole term given by 

Gmim 2 m3 2 


G 2 = 


16(mi + m 2 )a 2 (l — e 2 ) 3 / 2 


( 10 ) 


The semi-major axes and e 2 do not change at quadrupole 
order, so C 2 is also constant. We will henceforth refer to the 
dimensionless Hamiltonian, 


Hq = 


c 2 


( 11 ) 


2.3 Integrals of motion 

There are no dissipative forces in the problem, so the total 
energy, T-L, remains constant. Moreover, because no energy 
is transferred between the two binaries at quadrupole order, 
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Timescales of Kozai-Lidov oscillations 3 


each term of the Hamiltonian is conserved separately, so Hq 
remains constant as well. 

The total angular momentum is also conserved and may 
be expressed in the form of the geometrical relation, 


cos i = 


y^f2 S~l2 

^tot ~ ^1 ~ ^2 

2 G 1 G 2 


( 12 ) 


This relation is valid in the general case of three massive 
bodies. In the test particle limit the geometrical relation 
may be approximated by 


G t ot — G 2 + Gi cos i. 


(13) 


Since Gtot and G 2 are constant, we must have that Gi cosi 
is constant as well. Furthermore, Gi = L\j\, and Li is also 
constant, so this requires that jicosi be constant as well. 
We notate this constant of motion as 

0 = (jicosi) 2 . (14) 


In this form, the constant of motion is known as ‘Kozai’s 
integral’ and is equal to the square of the ^-component of 
the reduced angular momentum, j z (iHolman et al.l 119971) . 
Kozai’s integral implies that the component of angular mo¬ 
mentum perpendicular to the plane of the outer binary is 
constant. However, the test particle assumption is crucial to 
its derivation. In the general case of three massive bodies 
this component of angular momentum is not conserved, al¬ 
though it is possible to derive a generalized version which is 
conserved fe.g.. IWenll2003ll. The general i zed K ozai integral 
is due to the fact that, as Lidov fe Ziglinl (Il976l l showed, Hq 
is independent of g 2 , thereby implying that G 2 (and hence 
also e 2 ) must be constant. 

Because Hq only depends on ei, cosi, and gi, and there 
are two integrals of motion, Hq and 0, there is only one de¬ 
gree of freedom and so the system is integrable. Moreover, 
because these variables are all bounded, the motion is peri¬ 
odic (with the exception of a locus of stationary points of 
measure zero). The Hamiltonian to quadrupole order may 
be expressed as 

T-tq = \ [(5 - 3j1)(j1 - 30) - 15(1 -j1)(j1 ~ 0) cos 2<?i] 
Ji 

(15) 

in terms of j\ and 0. 


2.4 Equations of motion 


We are interested in the time evolution of the variables j i, 
cos i, and g\. Of these, only gi is a canonical variable so its 
time evolution follows directly from Hamilton’s equations: 


dg l 

dt 


BHq 


G 2 dHq 


dG i Li dj i 
Carrying out the differentiation of equation m we find 


(16) 


cf- = ^7 -T [5 ( 6 - A) (1 - cos 2gi) + 4 jf\ 


The variable j\ is related to a canonical variable, Gi 
constant, so we find its time evolution to be 


dj i 
dt 


1 OHq C 2 8Hq 


Li dgi Li dgi 


(17) 
by a 

(18) 


Again carrying out the differentiation of equation m we 
find 

(19) 

The time evolution of the inclination is complicated by 
the elimination of nodes from the Hamiltonian. Due to this 
procedure the time evolution of the inclination cannot be 
recovered from Hamilton’s equations directly. Instead, the 
inclination must be derived by calculating j\ and solving 
the geometrical relation given in equation m- 


2.5 Libration vs. rotation 


During a KL oscillation, the argument of periapsis of the in¬ 
ner binary may either rotate or librate. This is to say, g\ may 
sweep through the full range of angles from 0 to 27r (rotation) 
or it may be restricted to just a subset of them (libration). 
In the case of libration, the set of librating trajectories must 
librate about a fixed point of gi and ji . Inspection of equa¬ 
tion m reveals that j\ is stationary only when g i takes 
half- or whole-integer multiples of n (recall that 0 < ji). 
Now, inspection of equation m reveals that g i cannot be 
stationary at integer multiples of 7r. This implies that tra¬ 
jectories can only librate about half-integer multiples of 7r, 
so gi.flx = ±7t/2 and j 2 flx = ^50/3. 

To determine whether a particular system (i.e., a given 
H q and 0) librates or rotates we must see whether there 
exists a physical solution of equation m for j i when g% = 0. 
Setting gi = 0 in equation G3D and solving for ji we find 

ji =-A (10+ Hq +60). (20) 

The critical system on the boundary between libration and 
rotation will have a solution for ji exactly equal to unity 
and libration will occur if the only solution for j i exceeds 
unity. Defining the libration constant as 

Gkl = A- ( 2 - Hq - 60) , (21) 

we will have libration if Gkl < 0 and rota t ion if Gkl > 0. 
This constant was first presented in iLidovi dl962h and may 
be calculated equivalently by 

Gkl = e 2 ^1 - ^ sin 2 i sin 2 gi'j . (22) 

Note that the condition for rotation then becomes 


sin<?i ^ 


2 1 
5 sin* 


(23) 


Because Gkl naturally parameterizes a dynamical property 
of the triple, it is often convenient to work with it instead 
of Hq where possible. 


3 DERIVATION OF THE PERIOD OF KL 
OSCILLATIONS 

Because the Hamiltonian at quadrupole order is integrable, 
the period of KL oscillations, Ikl, may be determined ex¬ 
actly. The period may be written 

Ikl = dt = <j> ^j-dji. (24) 
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Solving equation m for cos 2gi and rewriting in terms of 
sin 2 gi, we have 


sin 2 gi = 


1 - 


'3 ji+jl ( 


H„ -90-5 +156 


15(l-Ji)0i-©) 


(25) 

Substituting equation (1251) into equation (1241) and substitut¬ 
ing the result into equation (1241) we find 


tj cl = 


30 C 2 


3 1 


1 - 


(!-Ji)0i-©) 

3ji + j\ (H q — 90 — 5) + 150 
15(1 — — 9) 


dji. (26) 


We note that this integral may be rewritten in terms of in¬ 
complete elliptic integrals of the first kind, but we do not do 
so here because it complicates the expression considerably. 

The integral in equation (1261) proceeds from the maxi¬ 
mum value of j i to the minimum value of ji and back again 
to the maximum value of ji , so we may instead integrate 
from j m in to jmax and multiply by two. Eliminating Hq in 
favor of Ckl by making use of equation (ED, and rearrang¬ 
ing, we have 


£kl = 


Li 

15 C 2 


1, 


1 - 


Jmin 

2 


(1 - Ji) 
1 0 


4 Ckl 
iF 5 1 -jl 


dji. (27) 


Eccentricity maxima (jmin) occur for gi = ±7t/2. Eccentric¬ 
ity minima (/ max ) also occur at <71 = ±7t/2 in the case of 
libration but occur at gi = 0 or n in the case of rotation. 
We may therefore solve for j m i n and j max by substituting 
the appropriate values of gi into equation m and solving 
for ji. We therefore have 


jmin = y i (C - VC 5- 600) 


(28) 


Jn 


(C + VC 5- 600), Ckl < 0 (29) 

(30) 


jmax — \/1 Ckl , Ckl 8 0. 

where we have defined 

() = 3 + 50 + 2C K l- (31) 

For convenience, we define the integral in equation ED to 
be /(Ckl, 0) such that 

154klC2 


/(Ckl, 0) = 


(32) 


Having calculated the limits of integration, we can now use 
equation ED to calculate the period of KL oscillations to 
quadrupole order in the test particle limit for any hierarchi¬ 
cal triple. 


4 A BRIEF SURVEY OF PARAMETER SPACE 

We now turn to a brief exploration of the range of values that 
the integral in equation ED may take. The overall timescale 


for KL oscillations is determined by the coefficient before the 
integral, which we present in more detail in Section T5. II The 
integral, however, depends on only two parameters describ¬ 
ing the triple: H q and 0, or equivalently, Ckl and 0. Thus, 
once the timescale of KL oscillations has been set, only two 
degrees of freedom remain. 

What values may Hq, Ckl, and 0 take? It is easy to 
see from equation m that 0 ^ 0 ^ 1 since both j i and 
cosi are bounded by 0 and 1. Moreover, it is clear from 
equation ED that 

— - < Ckl < 1 (33) 

since all the terms are bounded by 0 and 1. From the bounds 
on 0 and Ckl, we can conclude from equation ED that the 
bounds on Hq are 

- 10 ^ Hq < 20. (34) 

However, the limits on Hq and 0 are not independent. In 
the case of gi = 0 , the requirement that 0 ^ j\ implies that 

- 10 + 60 < Hq ^ 20. (35) 

This, in turn, translates to the requirement in Ckl that 


Ckl ^ 1 — 0. (36) 

In order for equation ED to have a solution, the square 
roots in equations (1281) . (1291) . and ESI) must exist. The ex¬ 
istence of the inner square root in equation ED requires 
that 

Ckl >-i (50 - 2V/150 + 3) . (37) 

This requirement is always satisfied in the case of rotation 
(Ckl > 0). In the case of libration (Ckl < 0), this require¬ 
ment may instead be written in terms of Ckl as 

0 ^ - (3 — 2%/—6Ckl — 2Ckl) , Ckl ^ 0. (38) 

In the case of libration, the square root in equation (1291) 
exists everywhere that the square root in equation (1281) does, 
so the existence of this square root adds no new constraints. 
In the case of rotation, the requirement that the square root 
in equation ED exist is satisfied by the same condition set 
in equation ED- 

Equation (1381) implies that there is a critical inclination, 
below which librating KL oscillations do not occur. Taking 
Ckl = 0 we recover the usual critical inclination of cosi ^ 
\J 3/5. Although we do not provide an explicit derivation, we 
note that the criterion in equation (1381) can also be arrived 
at by requiring that 



(39) 


when ji is at a maximum. In other words, KL oscillations 
occur when the minimum eccentricity is an unstable equi¬ 
librium. 

Knowing now the region of parameter space in which 
KL oscillations occur, we can numerically integrate the in¬ 
tegral in equation ED over the entire parameter space. The 
results of this procedure are presented in Figure [D Except 
for a narrow strip of parameter space centered around the 
boundary between rotation and libration (Ckl = 0 ) the 
integral only varies by a factor of a few. Near the rotation- 
libration boundary the integral diverges and KL oscillations 
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Ckl 


Figure 1 . Variation in the period of KL oscillations over all of 
parameter space. The contours (dotted lines) show different values 
of /(Ckl, ©) (he., the integral in equation I27j) . The period varies 
only by a factor of a few except very near the boundary between 
rotation and libration (dashed line), where it diverges. The gray 
regions indicate where KL oscillations are not possible. The large 
dot at Ckl = 0, © = 3/5 marks the largest value of © for which 
libration is possible. 


have arbitrarily large periods. Figure \l\ also indicates that 
the period of KL oscillations depends most strongly on Ckl 
and only weakly on 0. 


5 APPROXIMATIONS 


5.1 The timescale of KL oscillations 


So long as the integral in equation (ED is of order unity, 
the period of KL oscillations will be given by the coefficient 
before the integral to within an order of magnitude: 


Ikl — 


Lx 

15 CT 


(40) 


Substituting equations Q and m and noting that we are 
working in the test particle limit so m 2 —> 0, we have the 
timescale in terms of the semi-major axes, masses, and ec¬ 
centricities: 



This timescale may be expressed more elegantly in terms 
of the periods of the inner and outer orbits, Pi n and P ou t, 
respectively, by making use of Kepler’s law: 



This is the form of the KL period that typically appears in 
the literature, but with an additional mass term and numer¬ 
ical coefficient. The mass term implies that KL oscillations 
lengthen indefinitely as the tertiary approaches zero mass. 


In the case of a massive tertiary but a test particle primary 
and secondary (e.g., a WD-WD binary orbiting a SMBH), 
the period of KL oscillations approaches a constant value. 
Note that in the case of an equal mass primary and tertiary, 
neglecting the numerical coefficient will lead to an overesti¬ 
mate of the period of KL oscillations by a factor of nearly 
three. 


5.2 High inclination, low eccentricity triples 

In most cases of interest in astronomy, the inner binary of a 
hierarchical triple starts with a low to moderate eccentricity. 
Moreover, KL oscillations are strongest (and therefore most 
interesting) when the tertiary is at high inclination. It is 
therefore worth finding an approximation to tjcL in the high 
inclination, low initial eccentricity limit. In this limit, we 
have 0 —» 0 and Ckl —> 0 and equation (ZD may be solved 
exactly: 


/< c «-e)=4 X 

Jmin 

(43) 

We also have in this limit that j m i n 1 and 1 

so that 

Jmax 1 

/(Ckl, 0) — ~i= In (, 2 V 

4 V 6 \ 1 Jmax / 

(44) 

where 


Jmax — -L 1 ^ 

(45) 

for libration (Ckl < 0), and 


1 ^ KL 

Jmax — 1 — ^ • 

(46) 


for rotation (Ckl > 0). 

The dependence of /(Ckl, 0) on 0 is non-trivial to ap¬ 
proximate from first principles. After experimenting with 
several functional forms, we found that /(Ckl,0) varies 
most closely with (1 — 0). If there is a 0 dependence both 
inside and outside the logarithm, we then expect /(Ckl, 0) 
to take the form 

V0^) (i - er ' (47) 

where a, b, and c are fitting parameters. We fit numerical 
integrations of /(Ckl,0) to this form over the range 0 
0 0.25, —0.1 Ckl ^ 0.1 and find the remarkably good 

fit, 



We attempted to add several auxiliary parameters but found 
that they did not substantially improve the fit. 

The approximation provided in equation (1481) fits the 
true value of /(Ckl, 0) to within 2% over the range sampled, 
and over the vast majority of the range sampled the residuals 
are less than 0.3%. This is therefore an appropriate formula 
to use for triples in which the inner binary has an eccentricity 
ei < 0.3 and an inclination 60° < i < 120°. (Note that the 
triple need only have high inclination and low eccentricity 
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Ckl 


Figure 2. Residuals for the approximation in equation d48l) to the 
period of KL oscillations in the high inclination, low eccentricity 
limit. The approximation is correct to within 2% at all points 
in this range and typically does much better. This range of Ckl 
and © corresponds to triples in which the inner binary has an 
eccentricity ei < 0.3 and the inclination is i > 60°. 


at some point in the KL cycle for the approximation to be 
valid.) Contours of the residuals of equation (14811 are shown 
in Figured 


6 THE ECCENTRIC KL MECHANISM 

If the two masses of the inner binary are not equal and 
the outer orbit has non-zero eccentricity, the next term 
in the expansion of the Hamiltonian, the octupole order 
term, becomes dynamically significant. This term leads to 
changes to the orbital parameters of the outer orbit that 
are slow relative to individual KL oscillations. These long¬ 
term changes can cause the inner orbit to eventually pass 
through an inclination of 90°. During these orbital flips, the 
nearly perpendicular inclination leads to strong KL oscilla¬ 
tions which drive the inner binary to extremely large eccen¬ 
tricities. For this reason, the dynamical effect of the octupole 
term has been called the ‘ec centric KL mechanism’ (EKM) 
(e.g., iLithwick fe Nao3l201ll l. 

The introduction of the octupole term breaks the inte- 
grability of the Hamiltonian. Consequently, neither Ckl or 
0 remain constants of the motion. Furthermore, in the test 
particle limit at quadrupole order it is possible to eliminate 
the longitude of ascending node, fi, from the Hamiltonian. 
At octupole order either this parameter or g 2 necessarily 
enters into the equati ons of moti on . In this section we will 
follow the analysis of iKatz et al.l li’dll 1 1 ) and work in terms 
of the longitude of ascending node of the eccentricity vector, 
H e , defined such that e = e(sin i e cos fi e ,sin i e sin fi e , cos i e ), 
and e points toward periapsis of the inner binary. 

In the case of rotation (Ckl > 0), the parameters f2 e , 
Ckl, and 0 all change on a timescale which is long com¬ 


pared to individual KL cycles. It is therefore possible to 
assume that fi e , Ckl, and 0 are all approximately constant 
over individual KL oscillations and only examine the long¬ 
term changes to these parameters. In this approximation the 
system remains integrable with new integrals of motion. Due 
to the integrability of the system the variations in Ckl, 0, 
and are all strictly periodic. In this section we derive 
the period of these EKM oscillations. We note that there 
is a related octupole-order dynamical phenomenon in which 
nearly coplanar orbits at high eccentricity can undergo a 
flip by rolling over its major axis. An analysis of this phe¬ 
nomenon, including a derivation of the timescale, can be 
found in iLi et alJ ll2014blh 


6.1 Equations of motion and integrals of motion 


Since energy is conserved, the quadrupole order term of the 
Hamiltonian, T-Lq, is also conserved in the time-averaged be¬ 
havior of the system. This implies that the relationship be¬ 
tween Ckl, 0, and 7f q in equation (1211) remains valid and 
that the quantity 

<t> q = Ckl + i© (49) 


is a constant of motion. 

It is convenient to work with the parameter e oc t, which 
measures the relative size of the octupole order term of the 
Hamiltonian to the quadrupole order term. The parameter 
e 0 ct is conventionally defined as 


Coct 


e2 Q-i 

1 — e\ a ,2 


(50) 


Some authors have added a mass term (e.g., iNaoz et al.1 
1201331 1 to capture the fact that the octupole term is zero 
and EKM oscillations do not occur for an equal mass inner 
binary. However, because we are working exclusively in the 
test particle limi t we do not do so here. 

Following iKatz et al.l (faoill l. the long-term evolution in 
H e and 0 are given by 


dHe ( 6E{x)-3K(x) \ 

drr 4A'(a:) ) 


dO 157re oct \/@ sin fi e 

dr 64%/TO K (x) 


(4 — IICkl) \/6 + 4Ckl, (52) 


where K(x) and E(x) are complete elliptic functions of the 
first and second kind, respectively, 


, n \ _ 3(1 — Ckl) 
x(Ckl) = 


(53) 


3 + 2Ckl 

and the time coordinate has been scaled to the KL period 
during a flip: 

t 


(54) 


fKL,i=90° 

iKatz et~ahl 1 2011 1 also derive another integral of motion, 

X = F(Ckl) - £oct cos H e , (55) 

where the function F(Ckl) is defined to be 


F(Ckl) = 


32^3 


f 


K(r 7 ) - 2E( V ) 


n Jx(Ckl) — 21)y/2rj + 3 


dr/. (56) 


Although there are two integrals of motion, cj> q and x, 
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Timescales of Kozai-Lidov oscillations 7 


they are not sufficient to completely describe the dynamical 
behavior of the triple. This is because e oc t carries dynamical 
information as well, most importantly whether or not flips 
are possible. The dynamical significance of e oc t can be seen 
from the fact that e 0 ct enters into the definition of Thus, 
in the octupole case there are three independent parameters 
describing the system as opposed to the case of quadrupole- 
order KL oscillations in which there are only two. 


6.2 The period of EKM oscillations 

In the case of EKM oscillations it is easier to derive their 
period directly from the equations of motion rather than 
from action angle variables. We have from equation (|49jl 
that 


dC*KL 1 dO 

dr 2 dr ’ 


(57) 


so the period may be written 


tekm = 


dr 


dC KL 

Crl 


Substituting equation (1521) we find 


(58) 


term = 


128v70 K( x) 

157re oc t ^/2{(j> q — Crl) sin fi e 

X -- =dC < KL. 

(4 — 11 C'rl )'\/6 + 4 Crl 


(59) 


To write f! e in terms of Crl, we note that equation (ESI) 
implies that 


sin 


1 - 


X - ^(Ckl) 


Coct 


( 60 ) 


Substituting equation (1601) into equation (15911 and explicitly 
writing the limits of the integral yields 


6.3 Parameter space of the EKM 

As in the quadrupole case we first explore over what region 
of parameter space EKM oscillations with flips may occur. 
We then determine the variation in fEKM over this range of 
parameter space. Unfortunately, the parameter space cannot 
be mapped quite as easily as in the case of quadrupole KL 
oscillations because there are now three parameters describ¬ 
ing the system instead of two: cf) q , x, an d e 0 ct- As such, we 
explore parameter space for two choices of e oc t: £oct = 10 -3 
and e 0 ct = 10 -2 . Strong octupole-order effects occur in 
many triple systems with e oc t = 10 -2 , but these effects 
are much weaker for m ost triples when e oct = 10 -3 (e.g., 
iLithwick fc Nao3l201lll . 

To determine the boundaries of the parameter space of 
spin flips we first recall that 0 SC 0 ^ 1, and for rotation 
0 ^ Crl ^ 1 (which is the only case we are considering 
to octupole order). The occurrence of a spin flip is equiva¬ 
lent to having 0 = 0, and hence during a flip Crl = 4>q- 
Since cosfl e is bounded by ±1, we then have the following 
constraint: 


F((j> q ) - e oc t ^ X < F{<Pq) + «oct. (64) 

The parameter space can be divided into two regions based 
on the maximum of the function F((f> q ). This maximum can 
be found by solving K(x c nt) = 2F(x c nt ) for a; cr it, which 
yields a; cr it « 0.826, and then calculating 

Kent = w 0.112. (65) 

3 -I- 2 a: C rit 

Now, <j> q cannot be arbitrarily large because F(rf> q ) diverges 
at (j> q = 4/11. Thus we have 

K < YJ- (66) 

Since, for (f> q < (j> qt cnt , F{(j> q ) cannot be less than zero, this 
then implies a constraint on x : 

X^e oct {<Pq < 4>q ,crit ). (67) 


term 


256%/TO f cKL ’ m “ 


L 


K(x) 


157re oct JcKL.min V 2 (/> 9 - Crl) (4 - IICkl) 


1 - 


(x - F(Ckl))‘ 


(6 + 4Crl) 


dCyRL* ( 61 ) 


The upper limit of the integral can be deduced by noting 
that Crl is maximized when 0 is minimized and that 0 = 0 
during a flip. We therefore have 


©KL.max — ( 62 ) 

The lower limit is more subtle. It is clear from equation (1521) 
that 0 is maximized when sin 12 e = 0. This implies from 
equation EH> that 


E(C , KL,min) = X± eoct- (63) 

To decide whether to take the plus or minus sign, we must 
solve both for Crl and then take the value of Crl which 
is less than Crl, max- This equation can then be used to 
solve for Crl, min numerically. Together, equations 
and (16311 can be used to solve for the period of EKM oscil¬ 
lations exactly. 


Finally, the above relation implies that 


F{4>q, m in) = £oct- (68) 

Taken together, these relations bound the parameter 
space over which flips are possible. The resulting maps of 
parameter space for e oct = IIP 3 and e oc t = 10 -2 are shown 
in Figure [3] Because the parameter space over which flips 
occur is somewhat narrow and the dependence of term on 
( p q is fairly complicated we do not show contours as we 
did at quadrupole order in Fig. |T] Instead, we show term 
as a function of (j> q with the choice of x = F(tj> q ) and 
X = F((j) q ) ± t 0 ct/2 in Fig. [4] The timescale for EKM os¬ 
cillations depends most sensitively on rf) q . The timescale has 
two singularities: one at the maximum value of cj> q of 4/11, 
and another which is dependent on the choice of Xi but is 
near </ g , C rit ■ Except very close to these singularities, the pe¬ 
riod of EKM oscillations does not vary by more than a factor 
of a few. Thus, over a broad range of parameter space EKM 
oscillations have similar timescales. The existence of these 
singularities, however, does imply that Iekm has some de- 
pendence on, e.g. , the initial inclination as was found by 
iTevssandier et all (l2013l ~). 
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Figure 3. Parameter space where EKM oscillations with flips are possible for two choices of € oc t- We only explore the parameter space 
where individual KL cycles are rotating instead of librating (i.e., Ckl > 0), as librating cycles cannot be correctly analyzed using this 
technique of averaging over individual KL oscillations. At smaller values of e oc t the area of parameter space where rotating flips are 
possible shrinks about the line \ = E(^>). 




0q 0q 


Figure 4. The period of EKM oscillations with flips as a function of cf)q for three choices of %. The solid line is given by the choice 
X = F(cf)q ), the dashed line by \ = F(0q) — e 0 ct/2, and the dotted line by x — F(</>g) + e 0 ct/2. Except very near the two singularities, 
the period of EKM oscillations does not vary by more than a factor of a few. Over a broad range of parameter space EKM oscillations 
have similar timescales. 


6.4 The dependence on e oc t 

If the constants (f) q and x are held fixed and e 0 ct is varied, 
how does the period of EKM oscillations vary? Equation m 
exhibits a e oc t dependence in the coefficient before the in¬ 
tegral, so it is tempting to conclude that the timescale for 
EKM oscillations scales as e~^. This conclusion has been as¬ 
serted in several studies in the literature, but we show here 


that it is incorrect. The integral in equation m in fact 
exhibits a e~^f 2 dependence. 

To determine this dependence we first note that for 
EKM oscillations to occur, in general Ckl -C 1- This then 
implies that x is very close to unity, so we may write 
x = 1 — e, where 1. For values of x very close to unity, 
the complete elliptic integral of the first kind may be ap- 
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proximated 

K( 1 — e )~—line, (69) 

and the complete elliptic integral of the second kind is ap¬ 
proximated by _E(1 — e) ~ 1. We note that the coefficient 
in equation l|69|) is off by several tens of percent for realistic 
values of e, but the important feature of this approximation 
is that it carries the correct dependence on e. The function 
T’(C'kl) can then be approximated as 

F(CKL) ~ - J;V!rG ln(e ' )+2 ) *'• (70) 

Now, because we are integrating over a small range, the 
integral can then be approximated as 

j) + 2 )‘ < 7i > 

Furthermore, we note that 

£ ^ -Ckl (72) 

so we finally have 

f(Ckl) “ ~t&i/e Ckl Q ,n (tt) +2 ) ■ < 73) 

Let us now consider the lower limit of the integral in 
equation m- For simplicity, let us for the time being re¬ 
strict ourselves to the locus x = Fi&q) since here flips occur 
for arbitrarily small values of e 0 ct- We then have 

E(C K L,min) «* F(</> q ) ~ Coot . (74) 

Now, the approximation in equation m may be written 
more simply as E(Ckl) ~ kCn l, where fc is a parameter 
that has only a sub-linear dependence on Ckl. For small 
Ckl, then, the function F is nearly linear in Ckl - This then 
implies that for points on the locus we are considering 

- C-KL.min ~ (75) 

fc 

This then means that the width over which we are integrat¬ 
ing, ACkl is proportional to e 0 ct since 

ACkl = Ckl ,max Ckl ,min — (j> q — Ckl ,min ^ toct. (76) 

Let us now consider the various terms of the integrand 
of equation CD- We have already seen that because x is 
close to unity, K(x ) ~ 1ii(Ckl/3). This term is sublinear 
so we ignore it. The (4 — IICkl) term reduces to 4, and 
similarly the V6 + 4Ckl term reduces to \/6. The term 
\J2{4> q — Ckl) reduces by equation (l75l) to ~\/2e oc t. This 
leaves only the sin fi e term. Now, if <j> q — Ckl ~ £oct and F 
is approximately linear in this limit, it must be the case that 

F{4> q ) ~ Ckl = X — Ckl ~ £oct- (77) 




-to 


W 

-to 



Figure 5. The period of the EKM relative to the period of KL os¬ 
cillations as a function of e oc t calculated analytically using equa¬ 
tion ED (lines) and by integrating the secular equations of mo¬ 
tion (points). The timescale for the EKM is almost exactly pro- 
portional to e oct / and there is excellent agreement between the 
secular and analytic calculations. We show this relationship for 
an arbitrary choice of c()q = 0.015 and two choices of x — F{.<t>q) 
(black line and points), and x = F(<l>q) + 9 x 10 —4 (gray dotted 
line and points). For x — Fi&q) flips are possible at arbitrarily 
small values of e 0 ct, whereas for x = + 9 x 10 —4 flips are 

only possible for values of e 0 ct > 9 x 10 —4 . The relationship be¬ 
tween £ekm an d e oct becomes slightly shallower near this critical 
value of e 0 ct- Note that e oc t cannot exceed x■ Although flips oc¬ 
cur at larger values of e 0 ct, the evolution is no longer integrable 
because the inner binary switches between rotation and libration. 
In this regime the timescale for flips steepens as a function of e 0 ct, 
although there is no longer a simple relationship between the two 
because the evolution becomes essentially chaotic. 


Since the integral has a coefficient of e oc J., this then implies 
that the overall dependence of the period of the EKM is 


1 



(78) 


We demonstrate this dependence explicitly in Fig. [S] by 
numerically calculating the period using equation (16111 for 
fixed values of (f> q and x but over a range of e 0 ct- We have 
compared these values with the periods obtained by inte¬ 
grating the secular equations of motion directly and find 
excellent agreement. 

By combining this result with the numerical coefficient 
of equation (ED we find that 


Comparing this to equation we find that to lowest or¬ 
der, sin does not exhibit any dependence on e oc t. It is 
straightforward to verify this claim numerically. 

Putting these results together, we find that the only 
dependencies on e oc t in the integral come from the width of 
integration (which yields a dependence of e 0 ct), and from the 
term l/y/2((j> q — Ckl) (which yields a dependence of 


Iekm 


256VT0 

— -^=tKL,i=90° 

15tt yj £ oc t 


(79) 


During a single EKM cycle the inner binary will undergo two 
flips, so the flip timescale is half this value. The flip timescale 
can then be obtained by substituting for tKL,i= 90 ° i taking 
0 = 0 in equation (14811 . and noting from equation gi) that 
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Ckl = < j> q when 0 = 0, 


REFERENCES 


t n i i 


256 


97r 2 v /15e 0 ct 


1 + — 
m 3 


p 2 

out 


f, 2x3/2 

l-e 2 ) In 


9.42 \ 


V 'q 
(80) 

Over most of parameter space this expression is valid to 
within a factor of a few. For extremely large values of e 0 ct 
(e 0 ct ~ 0.1) our numerical experiments demonstrate that 
the dependence of /ekm on e oc t steepens and this expres¬ 
sion overpredicts the timescale for flips, but in this limit 
non-secular effects b ecome important, so the above analysis 




does not apply ( e .g., Antonmi__fcPgrgd 2012^ Katz_&Dond 


ap 

20121; ISetof[201E 


Antognini ct al 


lAntonini et al.l I20l4 Bode fe Wegd 12014 


2014ll . Moreover, the above analysis also re 


quires individual KL oscillations to be short relative to the 
EKM cycle. If this is not the case, then resonances between 
the quadrupole and octupole order te rms can ind uce chaotic 
variation of the orbital eccentricity (iLi et al.ll2014afl . Thus 
equation dH should be taken as a rough upper limit on 

tEKM- 


7 CONCLUSIONS 

Using action angle variables we have derived the period of 
KL oscillations at quadrupole order and in the test particle 
limit (equation 1271) . From the exact period we have derived 
the timescale for KL oscillations. We have explored the full 
range of parameter space over which KL oscillations are pos¬ 
sible and found that except very near the boundary between 
rotation and libration (|Ckl| <C 1) the period of KL oscil¬ 
lations does not vary by more than a factor of a few from 
the derived timescale (Fig.[TJ. By employing several approx¬ 
imations in the high-inclination, low initial eccentricity limit 
we have found a function that matches the true KL period 
to within 2% for triples for which ei < 0.3 and i > 60° 
(equation 148II . 

The strength of KL oscillations varies due to the oc¬ 
tupole term of the Hamiltonian. We average over individual 
KL cycles to calculate the period of EKM oscillations, and 
hence, the timescale for spin flips to occur. We map the pa¬ 
rameter space over which spin flips occur (Fig. [3]) and show 
that apart from near two singularities where spin flips do 
not occur, the timescale for EKM oscillations does not vary 
by more than a factor of a few (Fig. [4j. Finally, we show 
numerically and analytically that the dependence of e 0 ct on 
the timescale for EKM oscillations is (Fig. [5]) in con¬ 

trast to previous studies. We provide the EKM timescale in 
equation (® and the timescale for flips in equation ([501>. 
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